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Abstract 

Starting from the chemical master equation, we employ field theoretic techniques 
to derive Langevin-type equations that exactly describe the stochastic dynamics 
of the Frank chiral amplification model with spatial diffusion. The intrinsic multi- 
plicative noise properties are completely and rigorously derived by this procedure. 
We carry out numerical simulations in two spatial dimensions. When the inherent 
spatio-temporal fluctuations are properly included, then complete chiral amplifica- 
tion results from a purely racemic initial configuration. Phase separation can also 
arise in which the enantiomers coexist in spatially segregated domains separated by 
a sharp racemic interface or boundary. 



Introduction 



Questions about the origin of life have always held great fascination. One 
especially important problem is to explain chiral symmetry breaking in nature, 
why for example, it came to be that the nucleotide links of RNA and DNA 
incorporate exclusively dextro-rotary (D) ribose and D-deoxyribose while the 
enzymes involve only laevo-rotary (L) enantiomers of amino acids. Mirror 
symmetry is broken in the bioorganic world and life as we know it is invariably 
linked with homochirality. For recent reviews that survey existing hypotheses 
concerning this phenomenon, specific experimental realizations and additional 
references, see [1,2,3]. 

A simple chemical scheme to explain spontaneous generation of chiral asym- 
metry was proposed by Frank in 1953 [4]. In it original version, the Frank 
model leads to an amplified production of a chiral enantiomer from a racemic 
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mixture seeded with an initial bias. If however the initial state of the system is 
strictly racemic, it remains racemic for all time. The problem is to explain the 
origin of the initial chiral seed. The purpose of this Letter is to demonstrate 
that when the inherent spatio-temporal fluctuations are properly included, 
then complete chiral amplification results from a purely racemic initial con- 
figuration. Phase separation can also arise in which the enantiomers coexist 
in spatially segregated domains separated by a racemic interface or bound- 
ary. The internal fluctuations responsible for these two a-priori unpredictable 
outcomes are a natural consequence of the fact that chemical reactions are 
generally diffusion-limited, and take place in imperfectly mixed systems [5]. 

The general method we employ in this Letter is best appreciated by thoroughly 
working through a simple example based on an extension of the Frank model, 
which we now introduce. Let L and D denote a pair of enantiomers. Then 
the modified Frank model we will study is described by the following reaction 
scheme [6]: 

Autocatalytic production: 

ki ki 

L + A f L + L, D + Af D + D (1) 

K3 k 3 



Mutual destruction, or dimerization, in a second order reaction: 

L + D^P. (2) 



The ki denote the rate constants and we take the achiral substance A as a 
uniform constant background. The difference between this and the original 
Frank modejj] lies in the open-flow reactor nature of the process and the fact 
that the reaction (1) is allowed to be reversible (fc 3 > 0) [7]. The system is 
fed by an input of the achiral substrate A, whereas the output consists of the 
inactive product P, and the excess of the enantiomers. We further assume that 
each enantiomer diffuses with the same diffusion constant D and include this 
feature in the master equation description of this process. 



1 For the direct formation of chiral matter from achiral substrate, we would have 
had to include the steps A — > L(fc ) and A — > D(fco). According to [8], ko must be 
kept sufficiently small. Otherwise, these chirally unspecific reactions would proceed 
too rapidly thus generating large amounts of racemic matter that swamps the am- 
plification process driven by the autocatalytic (1) and mutual inhibition steps (2). 
To take this important observation into account, we simply set ko = from the 
outset. 
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2 Effective Bosonic Field Theory 



Our starting point for a systematic treatment of the above reaction scheme 
is an appropriate master equation. On a microscopic level, this comprises an 
exact description of the dynamics. From this equation, it is then straightfor- 
ward to derive an effective stochastic field theory. The salient steps are given 
in detail below. 

The chemical master equation for the kinetic scheme in Eqs.(l,2) is mapped 
to a "second-quantized" description following Doi [9]. Briefly, we introduce 
annihilation and creation operators and a J for L, bi and b\ for D at each 
lattice site z, obeying the commutation relations [a i5 aj] = 8 i j ) [6 i5 &j] = Sij. The 
vacuum state (corresponding to the configuration containing zero particles) 
satisfies a^|0) = 6j|0) = 0. Furthermore, a^|m^) = mi\mi — 1), a\\rrii) = |m^ + l) 
and a\ai\rrii) = m^rn^) (and similarly for the 6-sector). We then define the 
time-dependent state vector 

|*(t)>= £ P({m},{n},t)U(alr(bir\0), (3) 

{m},{n} i 



where P({m}, {n},t) is the probability distribution to find m^rii particles of 
type L,D, respectively, at each site z, and satisfies the master equation. Then 
by means of Eq.(3), the master equation can be written as a "Schrodinger" 
equation 

-T^W))- |*(t))=exp(- J fft)|*(0)), (4) 



where the lattice hamiltonian H is calculated to be 



H = % E (<4 - - «i) + f E (&! - &})(&« - &i) 

- E {ki((4) 2 ai - a\ ai ) + fci((6i) 2 6i - b\h) 

i 

+ fe 3 (a}a?-(a}) 2 a?) + *;3(6}6? 

+ fofai&i - aja^jfe.)} 5 (5) 

and (z,j) denotes the sum over all lattice sites z (with lattice spacing /) and 
their nearest neighbors j in a rf-dimensional space. 

Now take the continuum limit (/ — > 0) [10] in order to obtain the path integral 
representation of the underlying stochastic dynamics. Thus we can write the 
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time-evolution operator 

exp ( - Ht)= [ T>aVaVbT>bexp(-S[a, a, b, b}) 



(6) 



in terms of continuous fields a(x, £), a(x, i), 6(x, i), 6(x, i) with weight exp(— 5 [a, a, 6, 6]). 
Next, perform the shift a = 1 + a*, 6 = 1 + 6* on 5. This yields the shifted 
action for any space dimension d ) which is found to be: 

S[a, a*, 6, 6*] = J d d x J dt \a*(d t a - DV 2 a + k 2 ab - k r a + fc 3 a 2 ) 

+ 6* (d t b - DV 2 b + k 2 ab - k r b + fc 3 6 2 ) + a* 2 (fc 3 a 2 - ^a) 
+ 6* 2 (fc 3 6 2 - fcift) + fc 2 a*6*a6} . (7) 

The action S contains all the dynamics of the reaction scheme (1,2) including 
diffusion, and apart from taking the continuum limit, is exact. In particular, 
no assumptions regarding the form of the noise are required. At this stage, 
there are at least two independent ways to extract detailed information from 
S: by (i) employing field theoretic renormalization group techniques [11] or 
(ii) by exploiting the exact equivalence of this action to coupled Langevin 
equations. The former is suited for understanding chiral symmetry breaking as 
a problem in dynamic critical phenomena, the latter is ideal for understanding 
chiral amplification, evolution of enantiomeric excesses and net reactor yield 
as a problem in spatial pattern formation. 



3 The Exact Langevin Equations 

For the final step we use the gaussian transformation [12] which allows us to 
integrate exactly over the conjugate fields a*, 6*, appearing in the path integral 
/ VaVa*VbVb* e~ s ^ a ^ b ^ . This final step yields a product of delta-functional 
constraints which in turn, immediately imply the following set of exact coupled 
stochastic partial differential equations: 

d t a = DV 2 a + k\a — k 2 ab — k<$a 2 + /y a (x, t), (8) 
d t b = DV 2 b + fcxfe - k 2 ab - k 3 b 2 + %(x, t), (9) 

where the intrinsic multiplicative reaction noise is completely and exactly spec- 
ified as follows, namely (rj a ) = (775) = and 

(r/ a (x, t)7/ a (x', 0) = (fci - Mx, *))a(x, t) (S d (x - x')(f (* " 0> ( 10 ) 



4 



(rft(x, t) Vb (x.', t')) = (h - fc 3 6(x, t))6(x, t) 5 d (x - x')5(t - t'), (11) 
<?7 a (x, t) Vb (x', t')) = -^ 2 a(x, *)6(x, f) <5 d (x - x')5(t - t'). (12) 



In the absence of noise (mean field limit) the fields a(x, i), 6(x, £) correspond to 
the coarse-grained local densities of the L,D particles. With noise, these fields 
are generally complex, and do not represent the physical densities. However, 
the averages (a(x, £)), (6(x, t)) are real and do correspond to particle densities 
[13]- 

For convenience, we cast this system in terms of nondimensional couplings and 
fields. Introducing L and T for spatial and temporal dimensions, dimensional 
analysis implies that [a] = [b] = 1/L d , [D] = L 2 /T, [fa] = T~\ [k 2 ] = [k 3 ] = 
L d /T, [77] = l/TL d . So, define dimensionless fields a = {k 2 /ki)a ) b = (k 2 /ki)b, 
time r = kit, coordinates xj = ^JhjDxj and dimensionless noises fj a = 
{k 2 /k\)r] a ^ fjb = (^2/^)775. Then the equations Eqs.(8,9) can be written as 



d r a = V 2 a + a — ab — ga 2 + fj a . 
d T b = V 2 b + b-ab- gb 2 + fj bl 



(13) 
(14) 



and the correlation matrix of the dimensionless noises is given by 



B = e 



a(l — ga) —\ah 



\ -\~ab b(l-gb) 



(15) 



where g — |^ and e = f^^) is the spatial dimension-dependent noise 

amplitude. Note that in d = 2 this reduces to e = ^, indicating that the noise 
is controlled via a competition between the dimerization reaction and spatial 
diffusion. 

We apply the Cholesky decomposition B = MM T to extract the square root 
of the matrix in the form of a lower triangular matrix M [14]. We will use 
this decomposition to relate the noise to a new real white Gaussian noise £ = 
(£i, 6), (6(x,T)0(x , ,r / )) = ^•5(x-x , )5(t-t / ), such that (fj a ,fj b ) = r\ = Mf 
(and thus fj = f^M 7 ). Notice that by doing so the condition < fjfj >=< 
M^ ,T M T >= M < &f /T > M T = MM T = B is satisfied. This decomposition 
will allow us to separate the real and imaginary parts of the noise, a useful 
feature to have for setting up a numerical integration of the (complex- valued) 
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stochastic reaction-diffusion system [14]. We verify that 



f VWi o x 



M 



b 2 i 

( 



B11 / 



(16) 



(17) 



V 2 



pVa(l-^a) 1/2 (l-#a) 1/2 VWC 1 ~ fl^X 1 ~ 9^) ~ \^ 



V2, 



yields the desired decomposition of 5. 



4 Numerical results 



The non-trivial mean-field solutions of Eqs. (13,14), in the absence of fluctu- 
ations (deterministic case, e — 0), have been worked out in [6]: 

• (I) if g < 1: b — 0, a = ^, homochiral absorbing state solution, 

• (II) if g < 1: a = 0, 6 = ^, homochiral absorbing state solution, 

• (III) if(7>l:5 = 6= racemic active state solution, 

where <? = In solutions (I) and (II), one of the two enantiomers is amplified 
to its maximal value and the other one disappears, this corresponds to pure 
chiral amplification. On the contrary, in solution (III) the two enantiomers 
coexist with the same concentration, this corresponds to the racemic solution. 
This scenario is only valid in the mean-field limit. If the full problem is treated, 
then the noise due to spatial fluctuations and the diffusion of reactants must 
properly be taken into account. 

We will solve the full stochastic two-dimensional version of Eqs. (13,14) nu- 
merically, using reflecting boundary conditions and a finite difference scheme 
with At = 0.05, Ax = Ay = 0.51, and a grid of size L x L = 154 x 154 0. 

We first considered the evolution of a system with g = ^| = 0.643 < 1. Two 
representative results are shown in Figs. 1 and 2. The initial condition was 
set to a homogeneous, racemic situation (5 = 0.5, 60 = 0.5), and the noise 

2 In our numerical studies, the system size did not have any effect on the stationary 
solutions, only on the time required to reach equilibrium. 
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strength e = ^ = 10 -3 was fixed. A sequence of simulations were indepen- 
dently executed: every time the same mathematical problem was solved, the 
only difference being the random sequence of applied noises, and the result- 
ing evolution of each realization. In the two sets of rows, we plot the time 
evolution of the spatial distribution of the real part of a(x,^, r) (upper row) 
and r) (lower row). From left to right, the spatial distributions are 

shown at r = 0,83, 166, 1111,8333 and 140000, respectively. When displayed 
in color, red represents a concentration of 0, purple represents a concentration 
of roughly 0.5, dark blue roughly 0.6 and green is the maximal concentration, 
1.555 (= i). In both simulations, after the transient time (r > 100, i.e. after 
the first two frames of the sequences in Figs. 1 and 2), we found that even 
if the initial state was homogeneous and racemic (see purple squares in Fig- 
ures 1 and 2), the spatial fluctuations rapidly drive the system to a pattern 
with marked spatial compartamentalization, where, locally, the system is pure 
in one of the enantiomers GO. This means that we can find coexistence of lo- 
cal environments with the homochiral solutions (I) or (II). In the boundaries 
between regions of different chiralities there are racemic fronts of type (III), 
which in the mean-field limit would have only been expected for g > 1. Note 
the fully complementary nature between the spatial patterns formed of the 
a and 6-sequences during their evolution. Once the system has reached this 
state, the different regions compete and merge and the system evolves to one 
of the following final states: 

• A) A homochiral, homogenous solution where one of the enantiomers has 
entirely vanished, see Fig. 1. 

• B) A mixed state with strong spatial segregation of two immiscible phases, 
where locally in each phase, one of the homochiral solutions is realized, 
separated by planar racemic interface, see Fig. 2. 

Notice that when the solution is homochiral (solutions (I) and (II)) the noise 
vanishes since the correlation matrix, Eq. 15, vanishes as well. Thus this final 
state is absorbing. The racemic fronts are fluctuating and active. 

We verified that the complex solutions of Eqs.(13,14), once averaged over 
the fluctuations, do indeed yield real results, in accord with the theoretical 
expectations [13]. Next, we evaluate the time evolution of the spatially av- 

3 Thus from the third square onward from left to right in both Figures 1 and 2, we 
see the formation of red and green zones. Top row Figure 1, red and green signal 
the absence of and maximal concentration of a, respectively. Bottom row Figure 1: 
red and green signal absence of and maximal concentration of 6, respectively. Final 
state in this sequence is zero concentration a (red square top row) and maximal 
concentration for b (green square bottom row). In Figure 2, the final state (right- 
hand most squares in top and bottom rows) represents half the domain filled by 
maximal a and the complementary half square filled by maximal b concentrations, 
respectively, each half domain separated by a racemic boundary line. 
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a 




Fig. 1. Simulation 1. Starting from a homogeneous racemic condition, the local 
random fluctuations induced by imperfect mixing drive the system to a homochiral 
solution. Time runs from left to right. See the text for details. 



a 




Fig. 2. Simulation 2. The system is solved again with the same conditions as in Fig. 
1. The new random fluctuations induced by imperfect mixing lead the system to a 
different final state with spatial segregation. Time runs from left to right. 



eraged (over x, y) densities for the two cases shown above and for two other 
simulations representative of the inverse situation (i.e., where the system is 
dominated by the other enantiomer). We represent in Fig. 3 the time evolu- 
tion of the averaged enantiomeric excess ee(r) = < ? > ~ < ~ > . When the solution 

° v J <a>+<6> 

is spatially homogeneous with a at its maximal value, then ee — +1, for the 
inverse situation where b is amplified maximally, then ee = — 1, whereas for 
the intermediate situation with spatial compartmentalization < ee < 1 if the 
region with a is greater than the one with 6, and — 1 < ee < for the opposite 
situation. All these outcomes were found with roughly equal probability over 
the total of the 20 numerical simulations performed. 



As for the evolution of a system with g = ^| = 1.555 > 1 the system con- 
verges to the homogeneous solution a = b = 1/(1 + g) with small fluctuations 
about this value, i.e. the racemic state (III). All the results shown appear for 
a wide range of fluctuation intensities (e = 10 -2 — 10 -6 ). 
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Fig. 3. Time evolution of the spatially averaged enantiomeric excess 
ee(r) = of four representative simulations of the same study case, starting 

from racemic conditions (not resolvable in the figure). Due to the inherent symmetry 
of the reaction, either enantiomer a or b can dominate. The final result varies from 
one simulation to another. One can have total, homogeneous, chiral amplification 
(Simulations 1 and 4) or local chiral amplification with a racemic front (Simulations 
2 and 3). 

5 Conclusions 



i 1 3 

Si mutation 4 



Simulation 2 
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Internal particle density fluctuations automatically and inevitably lead to chi- 
ral symmetry breaking in an extension [7,6] of the Frank model. The inclusion 
of this noise is carried out through a first-principles technique which allows 
us to map the chemical master equation to the corresponding coarse-grained 
Langevin equations. This multiplicative noise is rigorously determined by this 
procedure and is not an ad-hoc addition to the mean field equations. A marked 
advantage of the fully stochastic model is that no initial chiral seed nor ex- 
ternal chiral field [15] are needed: full homochirality is achieved starting from 
a strictly racemic initial condition. In the limit as e — > 0, we would likely 
be sensitive to round-off errors, which also seem to induce chiral symmetry 
breaking, as reported in [16]. In two dimensions, there is a second class of 
solutions (see Simulation 2 in Fig. 2) for which the final state consists of two 
enantiomeric pure regions separated by a racemic interface or boundary. Re- 
garding the solutions displayed in Figs 1 and 2, it is amusing to point out 
that Frank, in conceiving of what may be expected to happen in imperfectly 
mixed systems, predicted the formation of separate "colonies" of the two kinds 
of enantiomers bounded by racemic surfaces, with chiral homogeneity main- 
tained within each such colony. He also stated that any initial curvature in 
these boundaries makes the boundary move towards the side from which it is 
concave, a feature which we have also observed in our simulations. The con- 
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sequence of this is that any enclave of the one species will eventually shrink, 
and the other survives; this is what we observe in simulations of type 1 (see 
Fig 1). An exception to this rule occurs "when the two seas are connected 
by a strait. A different species could survive in each sea, with the (racemic) 
boundary running stably from cape to cape" [4] . This is indeed what happens 
in simulations of type 2 (see Fig. 2) 

An entirely different stochastic version of the Frank model was presented in 
[17]. Chiral amplification is also achieved starting from a racemic initial state. 
There the noise is incorporated in an ad-hoc fashion by assuming fluctuating 
kinetic constants hi. 
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